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A numerical method is devised to solve a class of linear boundary-value problems for one-dimensional 
parabolic equations degenerate at the boundaries. Feller theory, which classifies the nature of the boundary 
points, is used to decide whether boundary conditions are needed to ensure uniqueness, and, if so, which 
ones they are. The algorithm is based on a suitable preconditioned implicit finite-difference scheme, grid, 
and treatment of the boundary data. Second-order accuracy, unconditional stability, and unconditional con- 
vergence of solutions of the finite-difference scheme to a constant as the time-step index tends to infinity 
are further properties of the method. Several examples, pertaining to financial mathematics, physics, and 
genetics, are presented for the purpose of illustration. © 2011 Wiley Periodicals, Inc. Numer Methods Partial 
Differential Eq 28: 807-833, 2012 
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I. INTRODUCTION 


The purpose of this work is to describe a numerical method for solving some classes of 
one-dimensional linear parabolic problems for equations, typically degenerate, of the form 


du 

9 7 


d 2 u 'du 

= a{x) — +b(x) — , 
dx z ox 


x € (r u r 2 ), t > 0, 


dv 

~dt 


d 2 d 

— (a(x) u) - —(b(x) v), 
dx z dx 


x G (r u r 2 ), t > 0, 


( 1 ) 
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where u = u(x,t), v — i >(x,t), — oo < n < x < r 2 < +oo. The second equation is the 
formal adjoint of the first one. They are also called the backward Kolmogorov and the forward 
Kolmogorov equation, respectively, and play a basic role in probability theory, because they pro- 
vide an analytical description of stochastic diffusion processes. The coefficient a(x) > 0 is the 
diffusion coefficient and b(x) is the drift. As is well known, “potential terms” like c(x)u, c(x)v, 
if present, can be removed from Eqs. (1) and (2), by setting for instance u = p(x)w, where p{x) 
is any solution of a related ordinary differential equation, which in general can only be done 
numerically. Therefore, we confine our attention only to equations without potential. 

The natural function spaces in which to pose partial differential equations (PDEs) like those of 
(1) and (2) are C°([/'i . r 2 [) and L 1 ((r i , ri)), respectively, as suggested by the theory of diffusion 
processes. In fact, Eq. (1) describes the time evolution of the moments, whereas Eq. (2) yields the 
time evolution of the probability density, of a given stochastic (diffusion) process. 

Throughout this article, we assume a'{x) and b(x) continuous but not necessarily bounded 
on (r l ,r 2 ). The diffusion coefficient, a (x ) , is assumed to be strictly positive on the open interval 
(ri, r 2 ), while it might vanish at one or both boundary points. The interval (r l5 r 2 ) itself may be 
unbounded. The PDEs (1), (2) are termed singular if at least one of the coefficients, a(x) and 
b(x), is unbounded at the boundary. We call them degenerate if the diffusion coefficient vanishes 
at one or both boundary points. In these cases, such equations have also been called (generalized) 
Feller equations [1], 

Singular or degenerate diffusion problems arise, for example, in mathematical genet- 
ics [2-6], in the theory of wave propagation in random media [7-11], and in quantitative 
finance [12]. 

The numerical treatment of initial-boundary value problems for such equations is nontrivial, 
being affected by a special pathology. In fact, in some cases, no boundary conditions (BCs) can 
be imposed, or, more precisely, a unique solution does exist (in some function space) without 
imposing any BC. It is therefore crucial to decide from the onset whether BCs are needed to 
ensure uniqueness in some appropriate function space, and, if so, which ones they can be. 

In 1952, Feller [13] (simultaneously with E. Hille [14]) realized that the uniqueness properties 
of solutions to problems for the two equations in (1) and (2) may be different. He also established 
that, in general, BCs cannot be prescribed independently of the behavior of the coefficients, a(x), 
b(x), near the boundary. Feller made a classification of the boundary points, r u r 2 , that can be 
used to decide whether a BC should be imposed on a given boundary point and, if so, which one 
it can be. 

Many works have appeared over the last 60 years concerning the analysis of elliptic and par- 
abolic PDEs degenerate at the boundary, see Ref. [15] and references therein, and very recently 
[16]. On the other hand, in recent years, the numerical treatment of linear (and nonlinear) degen- 
erate parabolic PDEs, in one or more space dimensions, has been considered in connection with 
the so-called Black-Scholes models, which started playing an important role in financial mathe- 
matics in 1973. In that context, the diffusion coefficient is positive inside the space domain, but 
vanishes or diverges to +oo at the boundary, see Ref. [17-22], just to mention few contributions. 
However, the problem of deciding whether BCs are required to ensure uniqueness of solutions 
(in some function class), and, in the affirmative, which ones, seems to be missing in such a broad 
literature. Typically, certain BCs were actually imposed which fact seems to confirm that they 
were determined in advance somehow (apparently without resorting to Feller’s theory). Usually, 
the aforementioned models possess a unique solution without imposing any BC on the boundary, 
but the correct behavior of such solution, clearly, was already known or determined somehow 
and then used in some numerical scheme. The main problem addressed in all these articles was, 
rather, that of implementing an efficient algorithm on a suitably truncated space domain, which 
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originally was unbounded, and then determining appropriate “artificial BCs.” These have been 
called “transparent BCs” when the solution in the so-obtained bounded domain coincides with 
the solution on the entire (unbounded) domain. It would be interesting to extend the results of 
this article to nonlinear singular Black-Scholes equations, because-compared with the classical 
linear models — they provide more accurate option values by taking into account more realis- 
tic assumptions [23]. The existing Feller’s theory, however, as it is, is only applicable to linear 
equations. 

Thus, it seems that the basic problem of deciding in advance whether any BCs could (or should) 
be imposed on the boundary points where the diffusion coefficient vanishes or goes to infinity 
was never considered in view of the ensuing numerical treatment. This observation applies also 
to other fields, besides that of financial mathematics. One can stress that this is a basic problem 
recalling that even for the Cauchy problem for the heat equation on the line, when the spatial 
domain is truncated to a bounded interval, according to Feller’s theory the nature of the boundary 
changes. 

In case that no BCs are required, numerical schemes that do not impose such data should exist. 
For this case, Feller theory also provides the limiting behavior of solutions or that of their related 
flux, and one can enforce this behavior numerically, instead of imposing any BC. 

In Section II, we derive a suitable numerical scheme, which is essentially an implicit Crank- 
Nicolson-type finite-difference scheme following a preliminary transformation and an appropriate 
boundary treatment. Preconditioning is also used, being mandatory at least in some cases. Feller 
classification of the boundary points is reviewed in Appendix A. In Section III, we give a number 
of examples, pertaining to a variety of fields, including option pricing in financial mathematics, 
wave propagation in random media, and genetics. Details needed to establish the nature of the 
various boundary points have been collected in Appendices B, C, and D. In a short summary at 
the end, the main points of this article are emphasized. 


II. THE NUMERICAL METHOD 

The numerical method developed in Section IIB below is based upon a preliminary transformation 
of Eq. (1) into a quasi-self-adjoint form (Section IIA). The method is then a Crank-Nicolson-type 
scheme applied to the transformed problem, using an appropriate grid and a suitable boundary 
treatment. We focus on those problems where the boundaries are reflecting (vanishing flux), as 
one or both boundaries are reflecting when there is a unique solution to the initial- value problem 
without imposing any BC (see the property mentioned at the end of Appendix A). In Section IIB, it 
is shown that the method is unconditionally stable in an Zri-norm, and that the numerical solution 
tends unconditionally to a constant, which approximates u 0 0 := lim,_ >+no u(x,t ) with an error 
independent of the time-step size. Numerical methods for quite general singular parabolic prob- 
lems have appeared previously in the literature (e.g. Ref. [24]) but apparently do not guarantee 
unconditional convergence to steady state with time-step independent error. 


A. A Preliminary Transformation 

The numerical method exploits the fact that Eq. (1) can always be written in a self-adjoint form 
having no potential (undifferentiated) term. Upon introducing the function 


s(x) 



x G (r u r 2 ). 


( 3 ) 
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(see (A5) in Appendix A), it becomes 


du 19/ 9m \ 

— — 5 a — I . 

3 1 s dx \ dx ) 


(4) 


This transformation, in other contexts, e.g., for the Schrodinger equation, is known as a “gauge 
change.” The subsequent change of independent variable 


-f 


£ = £(*):= / s(r) dr 


(5) 


is well defined because six) >0 for r t < x < r 2 and transforms (4) into the self-adjoint form 


du 9 / 9 du 
5 a 


d t 9§ V 3$. 

It will be assumed that the length L of the ^-domain is finite, 

r r 2 

L := §(r 2 ) - ^ (r, ) = / six) dx < oo, 

Jri 


(6) 


(7) 


that is, that s € L 1 ((?q, r 2 )). Admissible solutions shall first of all be bounded, so that in particular 




r 

/ u z d% < L sup it 1 , 

f(ri)<5 <?fr2) 


(B) 


and second of all be smooth enough that upon integrating (6) we may write 


1 d 

2 dt 

Here, / is the flux function 


-f(r 2 ) 


du 


2*L ^ = ^~L 


9m 


(9) 


f := s a — , 

9£ 


( 10 ) 


cf. (A6). What is sought is a solution existing for all times and satisfying the “reflecting barrier” 
BCs 


lim /= lim / = 0, t > 0. (11) 

In fact, we are interested in the cases when /q and r 2 ai'e either both natural boundaries, or one is 
natural and the other is an entrance point, in which case the initial- value problem for Eq. (1) has 
a unique solution in C°([/'i , r 2 1) without imposing any BC (Appendix A). The limiting behavior 
displayed in (1 1) is, however, known to be enjoyed by such solutions. Such solutions were stud- 
ied by Feller [3] for a special singular diffusion problem. Existence for our simple class of linear 
problems with time-independent coefficients can be studied by considering the spectral transform 
of the spatial operator, as was done, for example, by Morrison et al. [8]. General existence and 
uniqueness theorems for singular parabolic problems with time-dependent coefficients have been 
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established by Kohn and Nirenberg [25], Jamet [24], Carroll and Showalter [6], Friedman and 
Schuss [26], and Schuss [5,27]. 

The solution to our problem is unique in C°([r 1 ,r 2 ]) and tends to a constant as f — +oo. 
Uniqueness follows directly from (9), (1 1), and the boundedness of u, in the usual way: 


1 d 

2 dt 



u 2 d% 


r^n) 

= — s 2 a 

h(n) 



d$ < 0 , 


( 12 ) 


so the difference between any two solutions starting from identical initial data remains zero. In 
addition, the L 2 -norm of any solution is nonincreasing and hence tends to a constant as t — > +oo. 
The equality in (12) then implies that du/d% tends to zero, so the solution itself tends to a constant. 
This constant asymptotic solution, denoted by 


Moo := lim u(£,t), (13) 

t — +oo 

can easily be expressed in terms of the data. From (7) and (8), it follows that «(£, t) is integrable 
over the spatial domain, while (6) and (11) give 

d C Hiri) 

— I ud% = 0. (14) 

df 7f( n ) 

Therefore, 

r%(r 2 ) r^(r 2 ) 

Lu oa = lim / ud^= u 0 (£j)di;, 

l^+oo J m) 

where « 0 (£) := w(£,0), and we have simply 

] rS(r 2 ) 

Uoo = — / U 0 ($)d%, (15) 

the average of the initial data with respect to the measure d% = s(x)dx. Thus, the situation here is 
much the same as that for regular (i.e., uniformly parabolic) problems with vanishing Neumann 
BCs, and in fact our numerical method will apply equally well to such problems. 

We remark that the method also applies directly to forward Kolmogorov (i.e., Fokker-Planck) 
equations of the form 


dv 

~dt 


d 

dx 



d 

dx 


(b(x)v). 


(16) 


in case b(x) is linear in x, b(x) = fix + y, for then the substitution v := u exp {—fit} brings (16) 
into the form (1) with b := ^ — b. Linear drifts b(x) do arise in applications (e.g., Ref. [3]). On the 
other hand, for probabilistic applications it is often just as well to consider the backward equation 
as the forward equation (see Ref. [10] Section IV. A). The backward equation is somewhat simpler 
to treat numerically because in self-adjoint form it does not have a potential term. This is one of 
our motivations for taking as a starting point (1) instead of (2). 
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B. The Algorithm 

It is essential to our approach to cast the original equation ( 1 ) in the form (4) , which is not much of a 
restriction because (3) can be evaluated numerically if not analytically. As a matter of practicality, 
however, we do not require writing the problem in self-adjoint form: inverting the transformation 
£ = f (x) in (5) to obtain s 2 a in (6) as a function of £ is rarely possible to do analytically and can 
be tedious to implement numerically. 

By discretizing directly (4) instead of (1), it will be a simple matter to devise a numerical 
method that respects discrete versions of (12) and (15). From this, will follow easily two main 
properties of the method, namely unconditional stability, and unconditional convergence to steady 
state with an error independent of the time-step size. 

It is no loss of generality to assume that the length of the x -domain is finite, for if not, then let 
y = y(x) be any smooth mapping of (ty , r 2 ) onto a finite interval, y(r 2 ) — y(ty ) < oc. Such a 
mapping can be chosen to be explicitly invertible, in contrast to the mapping £ = £(x) discussed 
in Section IIA. Under the change of variable y = y(x), (4) becomes 


3 u 1 3 / _ 3 u\ 

— = — — I s a — I , 
3 1 s dy \ dy J 

with 



Hence, the form of (4) is preserved and assumption (7) is still met: 


py(r 2 ) pr 2 

/ s(y)dy= / s(x)dx = L<o o. 
Jy(ri) Jri 


(17) 


(18) 


(19) 


We have again the same problem, but now on a finite interval. This mapping, taking an originally 
unbounded space domain into a bounded domain, is effective, according to what was observed 
in Ref. [28], as the solutions sought are well behaved at infinity, as prescribed by Feller’s theory; 
see all numerical examples in Section III. 

Thus, we introduce a grid 


Xi := r x + (i - 0 Ax, /' = 1,2,...,/, Ax := (r 2 - u)/7, (20) 

I being a given positive integer, so that x\ = ty + Ax/2 and x t = r 2 — Ax/2. The numerical 
approximation to u(xj,t„) = u(xi,nAt), having set t n := nAt, n = 0, 1,2, . . ., will be denoted 
by u" . A discretization of (4) with second-order accuracy in time and space is then 


.v«(< +l - <) 

= n[p i+i - w; +l ) - A- ' « +1 - l ) + A+i ( u Ui ~ <) ~ A i « - <- 1)] (2D 
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for i = 2, 3 — 1, where 


1 At 


fi 


2 (Ax) 2 ’ 

Si := s( Xi ), @ i± i := s(x i± i) a(x i± i), 

( 1 !\ 

x ;± i := n + I i ± - - - I Ax, i = 2,3 1. 


( 22 ) 

(23) 

(24) 


Note that x 3 / 2 = r, + Ax and x/_i /2 = r 2 — Ax; no quantities are evaluated at the singular end 
points x = ri, x = r 2 . The BCs (11) may be written as 

3 u 3 u 

lim s a — = lim s a — = 0, (25) 

x-*r+ dx x ->r~ 3X 

so that with second-order accuracy we may write 

Sl {u1 +1 - u'{) = fl[@y 2 («2 +1 - < +1 ) + @3/2^2 - «?)], 

s/(«? +1 - < +1 ) = ti[-^ V2 (u n j +l - u n + 1) - /J 7 _ 1/2 « - (26) 

Equations (21) and (26) represent a system of I linear equations in I unknowns. Upon 
introducing the diagonal matrix 


S := diag(s!,s 2 , ...,s 2 ). 


(27) 


the symmetric tridiagonal matrix 



^ ^3/2 —@3/2 

... \ 


— @3/2 @ 3/2 + @5/2 

— @5/2 ■ ■ ■ 

B := 

@5/2 

@1- 3/2 + @1- 1/2 ~ @1-1/2 


V 

— @1-1/2 @1- 1/2 / 


and the vector 

u n := (i u\,u n 1 ,...,u n I ) T , 

this system may be written as 

S(n' ,+1 — u n ) = -fjL B(u n+1 + /<"), 


or 


(S + txB)u n+x = (S- t±B)u n . 


(28) 


(29) 


(30) 


(31) 


First, we want to verify that the scheme is well defined, i.e., that the symmetric matrix S + n B 
in (31) is invertible. By direct calculation, one has 

/-i 

yT B v = ^ @ i+ i (u i+ i - v t ) 2 (32) 

/=1 
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for any vector 


v = (v u v 2 , ...,n/) r . (33) 

Thus, the symmetric matrix B is positive semidefinite, and its sole null vector is 

e:={\,\,...,\) T . (34) 

As B is positive semidefinite and S is clearly positive definite, it follows that S + /r B is positive 
definite and hence invertible. 

Unconditional stability is immediate, for by left multiplying (30) by ( u n+l + u") T and using 
(32) one has 

i i-\ 

E4« +1 ) 2 - «) 2 ] ^ -/*E A- + > [« I 1 + «7 + i) - {O 1 + Of < 0, (35) 


i = 1 


i = 1 


which is a discrete version of (12). In fact, this means that the £ 2 -norm of the vector u" is 
nonincreasing with n, as the L 2 -norm of the solution u is nonincreasing with time. In particular, 


IX « +1 ) 2 < 


(36) 


i = 1 


i= 1 


The scheme also conserves “mass,” analogously with (14), as follows. As all columns sums of 
the matrix B vanish, one has for any vector v that 


= °- 


(37) 


1=1 


From (31) it then follows that 


Si u’l +i — Si u", for all n > 0, 


(38) 


i=i 


i=l 


the discrete analog of (14). Finally, observe from (36) that the sequence 


(Q 2 | , n = 0, 1,2, . . ., 


(39) 


1 = 1 


is monotone nonincreasing. It is also bounded from below, by zero, and therefore has a finite 
nonnegative limit as n -> oo. From (35), we then have 


7-1 


iim yft + ,/2[« + + ,' + Oi) - (o 1 + of = o, 

n — ^oc : ■ - 


(40) 


1 = 1 


and therefore 


«+/ + < +1 ) - (u'; +l + <) ^0 as n 


oo, 


(41) 
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for each i = 1,2,...,/ — 1. That is, 

u n+1 + u" — > 2 u°° e as n — > oo (42) 

for some number u°°, where the vector e was defined in (34). As Be = 0 and S is positive definite, 
it follows from (42) and (30) that 


,«+ 1 


0 as n — ► oo. 


From (42) and (43), we obtain 


as n 


oo, 


(43) 


(44) 


the discrete analog of (15). Thus, as n — > oo, the numerical solution u" tends unconditionally to 
a constant, u°°. To evaluate u°°, write (38) as 


i i 

Y Si w” = Yj s i Ll< ! f° r a " n >0. 

i = 1 i=l 


Letting n — > oo and using (44) yields 


u 


OO 



(45) 


(46) 


the discrete analog of (15). Clearly, u°° is independent of At, and therefore so is the error «°° — u^. 

In the examples worked out in Section III below, where lateral conditions of vanishing flux 
have been implemented, we also solved numerically the PDEs satisfied by the flux v(x,t) := 
ci(x) s(x) u x , in order to check the monotonic behavior of the solution u, which was apparent. In 
fact, the function a(x) s(x) being strictly positive on the entire domain, the sign of the flux is the 
same as that of u x . Starting from the quasi-self-adjoint form (4), i.e., s(x) u, = v x , we obtain, 
differentiating with respect to x, 


dv ( 1 dv\ 

3 1 \s dx J 


(47) 


The initial condition is v(x,0) = a(x)s(x)u x (x, 0) = a(x)s(x)u' 0 (x), while the BCs will be 
v(rj,t ) = 0 for j = 1, 2. These are discretized as 

Sl (u n + 1 - <) = ii[Ps (wf 1 - m" +1 ) - fix < +1 + pi (u n 2 - ill ) - Piu'l] 

at the first grid point, j — 1 , and 


S] (u 


n + 1 n 


V = li[-P I+l u'} +1 -p I i{u* I 


n + 1 __ n + 1 
U I - 1 


) i( M " M "-i)] 


at the last grid point, j = / . The terms ft and P l+ 1 are evaluated at the end points and hence 


they might be singular. Expanding to first-order s . i ~ Sj + ^s'., we obtain the estimate 


3+2 


2 "3’ 


^±i a j± I ( 5 3 ± 


7/ Ax v' 

= <3 . , i ± a.,i — . 

3±2 3± 2 2 5 ; 
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The term sj with j — 1 or j = I is not singular, while the diffusion coefficient vanishes at the 
end points. Therefore, the ratio ft ± 1 /s n for j = 1 or j = 7, also vanishes. The discretization of 
the BCs thus becomes 

W" +I (si + lift) 3) — [lft 3 U" +l = — 11ft 3) + iift^u'ft 

u n , + 1 (s, + \xftj_ 1) - = u"(s, - lift,_i) + !ift I _iu n ,_ l , 

which are similar to those of vanishing flux. 

In the PDE (59) in Section IIIC below, the BC at p = 1 is given by u(\,t) = 1, which is 
discretized at the 7 th grid point as 


s, («';+' -u n j) 


= p[fti+ i( J 

-uT l ) 


— u 

1 ) + 


Hence, the matrix B becomes 






-ft. 

0 


0 


-h 


■H 

0 


B = 

0 




0 








l 0 


0 




-u n i- 1 )]- 


\ 




m 


(49) 


and the system to solve becomes 


(S + 11 B ) U n+1 = (S - 11 B) U" + C, 


(50) 


where C = (0, . . . , 0, 2[ift, 1 ) 7 is the vector of the known terms, all zero except that evaluated 
at the 7th grid point. 

Typically, the implicit finite-difference method developed above needs to be preconditioned. 
In fact, because of the degeneracy, one or both the first rows of the matrix in (28) are quasi singular 
when Ax is small and vanish as Ax — > 0. Therefore, ill-conditioned matrices would occur for 
the most trivial reason, namely that their determinant is very small. Resorting to the simplest 
preconditioning strategies, i.e., either using the diagonal as a preconditioner or using the incom- 
plete Choleski factorization, is shown to circumvent such a difficulty. For instance, the effect of 
a diagonal preconditioner on tridiagonal matrices is clear. In fact, assume that a (/g ) = 0. Then, 
premultiplying the tridiagonal matrix (7)y) ; j =1 m by P ', where P := diag(Tn, . . . , T mm ), takes 
the quasi singular terms T n and T l2 into 1 and T x 7^12 , respectively, and the latter are both of order 
1, as Tn and T ]2 are both of order Ax. On the other hand, the tridiagonal structure is preserved. 

We have applied both preconditioners, the first to the system for the financial mathematics 
model and the second to the model of wave propagation in random media. 


III. APPLICATIONS 

In this section, we present a few examples, pertaining to a variety of areas. 
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A. Financial Mathematics 

In quantitative finance, the time evolution of the contract price, u(S, 1 ), of a “European option,” 
is governed by the linear parabolic one-dimensional equation 

u, = S 2 u S s — ( aS — b)u s , 0 < 5 < +oo, 0 < t <T, (51) 

where a > 0 and b > 0 are two given constants [12]. Such equation is very close to the so-called 
Black-Scholes equation, the potential term being missing here. Incidentally, a European option-as 
opposed to an American option — is an option that can be exercised only at the end of its maturity. 
The corresponding PDE problem for pricing European options is simpler, being set on a fixed 
domain, whereas the American options involve a free-boundary problem. 

The “underlying” asset price 5, here an independent space-like variable, is interpreted, in the 
underlying probabilistic model, as a stochastic mean reversion process, S(t,a> ) (to denoting the 
chance variable), obeying the stochastic differential equation dS = a(L — 5) dt + a dW . Here, 
L represents the long time average, i.e., the price at which 5 will tend to stabilize after long times, 
a and a are strictly positive constants, representing the mean reversion speed and the process 
volatility, respectively, and W is the one-dimensional standard Brownian motion. 

Equation (51) is of the type of Eq. (1) and is degenerate at 5 = 0 and singular at 5 = +oo. 
The boundary point 5 = 0 turns out to be an entrance point, whereas 5 = +oc is a natural point 
(see Appendix B). 

Feller theory then guarantees that the initial-value problem for the PDE in (51) has a unique 
solution in C°([0, +oo]) without imposing any lateral condition. Note that such theory, yielding 
uniqueness in the closed interval [0, +oo], implies boundedness on the halfline. 

Feller theory also provides some lateral conditions, which are properties of the aforementioned 
unique solution. These are those of vanishing flux, which, in case of Eq. (5 1), is the function 



Incidentally, note that 

lim f(S,t)= lim f(S,t) = 0 (53) 

S-J-0+ S— >+o o 

whenever u s is bounded (or diverges suitably slowly) near the boundaries. In this case, we would 
recover directly the properties of vanishing flux as prescribed by Feller theory. Indeed, we can 
prove that u s is bounded near the boundary 5 = 0. If u solves the problem above, differentiating 
both sides of the PDE in (51) as well as the initial value with respect to 5, it appears that the 
function v := u s solves the problem v, = S 2 v S s + [b — (a — 2)5]u s — av, v(S, 0) = u' 0 (S). The 
potential term, —av, can be promptly removed setting v(S, 1 ) := w(S, t) e~ at . Thus, w solves the 
problem w, = S 2 w S s + [b — (a — 2)S\w s , w(S, 0) = u' 0 (S). We can now apply Feller theory to 
such problem, first classifying the boundary points. Note that only the drift term differs from the 
PDE satisfied by u, being the same with a — 2 in place of a. As a typical value of a is greater than 4 
[12, Table 1, p. 340], we leave to the reader the analysis for 0 < a < 4 and consider only the case 
a > 4. With a — 2 > 0, all results established for the problem satisfied by a also hold for w , and 
we can conclude that the boundaries have the same nature as in the problem for u. Consequently, 
u s (S,t) is unique and continuous up to and including the point S = 0. In particular, u s is bounded 
in the right neighborhood of S = 0 and the BC of vanishing flux, S~ a e~su s = 0 is correct, as 
S- a e~-s =0 for 5 = 0. 


Numerical Methods for Partial Differential Equations DOI 10.1002/num 


818 


CACIO, COHN, AND SPIGLER 


European option value 



FIG. 1. European option price at the initial time, t = 0, and at the maturity time, t = T = 6 years. This 
is shown for asset prices S e [0, 100]. [Color figure can be viewed in the online issue, which is available at 
wileyonlinelibrary.com.] 

We can go further, and obtain an estimate for u ss , showing that it is bounded near 5 = 0. Now 
z := u S s satisfies z, = S 2 zss + [b — (a — 4)5]zj + 2(1 — a)z, z(S , 0) = Wq( 5). Note the favorable 
circumstance that the drift in (51) is linear implies that a term proportional to u s does not appear 
in the latter PDE. Again, the potential term can be removed by setting w(S, t) := z(S, t ) e~ 2(a ~ 1)l , 
obtaining w, = S 2 w ss + [b — (a — 4)5 ]»,><,, w(S, 0) = m ( "( 5). Finally, applying Feller theory once 
more, the boundary points turn out to be again of the same kind ( a >4), and hence w is unique 
in C°([0, +oo) x [0, T |). In particular, uss is continuous up to the boundary 5 = 0 and also 
bounded there. Therefore, the diffusion term S 2 u S s tends to zero as 5 — 0 + . We could say that at 
5 = 0, the PDE ceases to be parabolic and becomes formally hyperbolic, and the sign of the drift 
at 5 = 0 shows that the characteristics point outwards from the domain, justifying why there is 
no need of a lateral condition to make a solution unique. In Fig. 1, the solution w(5, t) to problem 
(51), i.e., the contract value of a European option, is shown at the maturity time T = 6 years, in 
the range of the asset price 5 e [0, 100]. This has been computed by the Crank-Nicolson-type 
finite-difference method developed in Section II, with I = 1000. Note that the solution decreases 
to zero in the right neighborhood of 5 = 0, while the option price stabilizes for 5 large. This 
phenomenon is compatible with the financial model being analyzed, because in a mean-reverting 
process, the asset price, 5, stabilizes over long times around a mean value, and consequently the 
option value tends to a constant value as well. 

It is worth noting that, in Ref. [12], the same problem was tackled by a semispectral method. 
The authors there transformed the original PDE in (51) into 

U, = x 2 U xx + x (a + 2 — x) U x , 0 < x < +oo, 0 < t < T, (54) 

setting x = b/S, and thus exchanging the boundary points 0 and +oo with each other. The PDE 
in (54) can be written in the form 


U, = 


1 

— -[x w(x) U x L, 
w(x) 


( 55 ) 
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Evolution in time of option value 



Asset price S do Time t 


FIG. 2. Time evolution of the option value for S e [1,6], with a maturity time T = 6 years. Here, the values 
have been computed by the semispectral method developed in Ref. [12], using 100 Laguerre polynomials. 
[Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 


where w(x ) := x a e~* is the weight function of the generalized Laguerre polynomials. This sug- 
gests writing the PDE (55) in abstract form and expanding the solution in a series of generalized 
Laguerre polynomials (that will be suitably truncated later). In Fig. 2, we have plotted the time 
evolution of the option value for S e [1,6], with a maturity time T = 6 years. The values have 
been computed by the semispectral method developed in Ref. [12] (compare with Ref. [12], Fig. 1, 
p. 341), using 100 Faguerre polynomials. In this case too, the solution decreases near x = 0 and 
stabilizes for large S, after long times. For values of S close to the origin, however, the method 
yields rather unclear results, because the solution oscillates unboundedly, see Fig. 3 (cf. Ref. [12], 
Fig. 3, p. 342). This behavior represents a phenomenon incompatible with Feller theory, which 
predicts — in this case — uniqueness and continuity of solutions in the closed interval [0, +oo]. 
The picture is meaningless also from the financial point of view: when the asset price tends to 
zero, the contract value of the option is also expected to decay to zero (rather than oscillate indef- 
initely between — oo and +oo). Note that here time, measured in years, goes backwards, from 
the maturity time, T, to the initial time, see Ref. [12]. Applying directly Feller classification to 
the PDE in (54) shows, by standard calculations, that the boundary points have the same nature 
as for the PDE in (5 1). It follows that, in contrast to what is shown in Fig. 3, the unique solution, 
U (x, t), is indeed continuous up to and including the boundaries. 

B. Wave Propagation in a Lossy Random Medium 

In the theory of wave propagation in lossy random media, the evolution of the moments of the 
complex- valued reflection coefficient is governed by the equation 


u, = -(1 - p Yu pp + 


1 , ,1 

t-Q-p )“ - 

4 p 2 


0<p<l,f>0, 


(56) 
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Oscillations In option values for small values of S 



Asset price S 


FIG. 3. Option price at 6 years from the maturity time, evaluated by the semispectral method of Ref. [12] 
using 100 Laguerre polynomials, for asset values in the range [0, 1], as shown in Ref. [12]. [Color figure can 
be viewed in the online issue, which is available at wileyonlinelibrary.com.] 


with u(p, 0) = p 2m ,m = 1,2, . . ., where a > 0 is the lossy parameter [9]. Here, t does not 
represent time but the width of the slab occupied by the random medium. The quantity u( 0, t ) 
then yields the m-th moment of the reflected power. 

Equation (56) is of the type (1) and is degenerate at p = 1 and singular at p = 0. The singularity 
can be removed by setting x := p 2 , obtaining 

z t = x (1 — x) 2 z xx + [(l-x) 2 -o'i]t I , 0 < jc < 1, t > 0, (57) 

with z{x, 0) = x'",m = 1,2,..., which is degenerate at both boundary points. The nature of such 
points is the same in both cases, and x = 0 (as well as p = 0) is a natural boundary, while x = 1 
(as well as p — 1) is an entrance point, see Appendix C. 

Also in this example, the solution is unique in the class C°([0, 1]), and such solution has the 
property of having the flux vanish at the boundaries. This is the condition that will be imple- 
mented numerically. In Fig. 4, we show the time evolution of the moments of the complex- valued 
reflection coefficient for a wave propagating in a lossy one-dimensional random medium, with 
the parameters m — 1, a = 0.5, and I = 1000. The solution decays in time and tends to attain 
a stationary profile, as predicted by the theory. Hence, for large 1 there is a stationary solution 
corresponding to the ordinary differential equation u, = 0. Note the detailed monotonic behavior 
of u(x, t) in the left neighborhood of x = 1. This has been confirmed from the sign of u x , which 
has been obtained by solving numerically the parabolic equation satisfied by the flux, see Fig. 6. 
See also Fig. 5, where the parameters are m = 2, a = 0.5, and I = 1000. 

In Fig. 7, the mean reflected power, with m = 1, a = 0.5, and I = 1000, is shown. The picture 
shows that, in a sufficiently thick slab filled in with a one-dimensional random medium, the mean 
reflected power attains a constant value. From the physical standpoint, in fact, when the slab is 
sufficiently thick ( t is large enough), a constant value of power is reflected back, i.e., the mean 
reflected power, u( 0, t), is constant. The same occurs to the mean dissipated power. 
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Evolution of the moments of the reflection coefficient (with dissipation) 



FIG. 4. Time evolution of the moments of the complex- valued reflection coefficient for a wave propagating 
in a lossy one-dimensional random medium, with m = 1, a = 0.5. [Color figure can be viewed in the online 
issue, which is available at wileyonlinelibrary.com.] 


C. Wave Propagation in a Lossless Random Medium 

When the random medium is ideally lossless, the same quantity, u, described in the previous case, 
obeys instead the evolution equation [7] 


u, = (1 - P 2 ) 2 u pp + 


i(l -P 2 ) 2 
P 


Up, 


0 < p < 1, t > 0. 


(58) 


Evolution of the moments of the reflection coefficient (with dissipation) 



FIG. 5. Time evolution of the moments of the complex-valued reflection coefficient in the lossy case, with 
m — 2, a = 0.5. Also in this case, the solutions stabilize over long times, as suggested by the theory. [Color 
figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 
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Flux in wave propagation with dissipation 



FIG. 6. Time evolution of the flux for the problem of waves propagating through a one-dimensional lossy 
random medium, with a = 0.5, m = 1 . This flux turns out to be positive at each point and each time, hence 
confirming that u grows with x, see Fig. 4. [Color figure can be viewed in the online issue, which is available 
at wileyonlinelibrary.com.] 


This equation coincides (rescaling “time”) with Eq. (56) with a — 0 and can be rewritten as 


du 

~dt 


(1 -P 2 Y 3 


du 


p — 
dp |_ 9p 


(59) 


Mean power reflection (with dissipation) 



Time t 


FIG. 7. Plot of the mean reflected power versus the slab thickness, with m = 1, for a wave propagating 
in a one-dimensional lossy random medium, with a = 0.5. [Color figure can be viewed in the online issue, 
which is available at wileyonlinelibrary.com.] 
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and setting 


x = 


1 + p 2 

1 -p 2 ’ 



in the conservative form 


3 u 1 3 f , 3 u 

— = (x 2 — 1 ) — 

3/ 4 dx dx 


1 < x < +oo. 


(60) 


An explicit formula for its solution can be obtained in terms of the so-called Mehler transform, 
see Ref. [7], but this formula seems to be very hard to use in practice, even through its numerical 
evaluation, being affected by a number of singularities. 

Equation (58) is degenerate at p = 1 and singular at p = 0; p = 1 is a natural boundary, while 
p = 0 is an entrance point, see Appendix C. We can prove that when p — > 1 ~,u p and u pp remain 
bounded, hence Eq. (58) reduces to u, = 0 there. As p = 1 is a natural boundary, we can impose 
there the condition of vanishing flux, 


lim /(p) = lim pu p = 0, 

p-> i _ p-> i~ 

hence w p | p=1 = 0. Thus, u p is bounded near p = 1 and we can impose the lateral condition 
u(l,t) — const. Integrating directly the equation u, = 0 and using the initial value «(p, 0) = p 2m 
(to obtain the /nth moment of the reflection coefficient), we obtain 


u(l,t) = u( 1,0) = u(p, 0)| p=1 = 1. (61) 

Note that the coefficients of both equations, (56) and (58), are “symmetric” in p around p = 0, 
hence, one could expect that their solutions have the property u p \ p=0 = 0. Such a property comes 
from the nature of the polar coordinate p. In some cases, one can obtain a symmetric problem 
taking p = 1 into the origin and proceeding similarly. This was done in Ref. [29,30] and [1 1]. 

In the lossy case [Eq. (57)], Feller theory then guarantees uniqueness of solutions in the class of 
continuous functions C°([0, 1] x (0, T |). Moreover, the boundary point x = 1, being an entrance 
point, is a reflecting boundary, so we can impose at x — 1, in the numerical scheme, the condi- 
tion of vanishing flux. At the natural boundary x = 0 we can also prescribe a reflecting barrier 
behavior, see Appendix C. The same can be said for the formulation in (56): the boundary point 
p = 0 is again an entrance point, thus we can impose there the condition of vanishing flux; the 
natural boundary p = 1 is also a reflecting barrier. 

The flux / for Eq. (58) is 


fix, t ) := — exp < —a 
x 0 


1 — X 1 — Xq 


U x . 


and lim t _ >0 + f(x, t) = lim,^ r fix, t) = 0, as predicted by Feller theory at natural and entrance 
boundaries. Therefore, \2*Jx u x \ x=0 = 0, i.c., u x (x, r) is bounded or diverges more slowly than the 
vanishing rate of »Jx as x — > 0 + . Therefore, the flux also vanishes as x —>■ 0 + and this condition 
can be imposed at rj = 0. 

In Fig. 8, we show the time evolution of the moments of the complex-valued reflection coef- 
ficient, with the initial data u{x, 0) = x 2 and the parameters m = 1 and I = 1500, in case of a 
wave propagating through a one-dimensional lossless random medium. The picture shows that 
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Evolution of the moments of the reflection coefficient (without dissipation) 



FIG. 8. Time evolution of the moments of the complex-valued reflection coefficient with the initial data 
u(x. 0) = x 2 , with m = 1, for a wave propagating in a one-dimensional lossless random medium. [Color 
figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 

the solution attains a stationary profile after long times. Hence, also in this case, we obtain a 
stationary value for u (p, t), for large t, correspondingly to the solution of the ordinary differential 
equation u, = 0. 

In Fig. 9, the mean reflected power is depicted as a function of the slab width, with m — 1, 
/ = 1500, for a wave propagating in a one-dimensional lossless random medium in the slab. The 

Mean power reflection (without dissipation) 



Time t 


FIG. 9. Plot of the mean reflected power versus the slab thickness, with m = 1, for a wave propagating 
in a one-dimensional lossless random medium. [Color figure can be viewed in the online issue, which is 
available at wileyonlinelibrary.com.] 
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picture shows that the mean reflected power stabilizes to a constant value. In fact, when the slab 
filled in by the random medium is sufficiently thick (i.e., for t sufficiently large), a constant value 
of the power is reflected, i.e., the mean reflected power, u( 0, t), tends to a constant. 

D. Genetics 

The time evolution of large populations of independent individuals, whose reproduction rate is 
independent of the size of the population, is described by the linear Fokker-Planck equation 

v t = (axv) xx — [(bx + c) v] x , 0 < x < +oo, t > 0, (62) 

where a > 0, b, and c are given constants [3]. 

Unlike those of the preceding examples, Eq. (62) is of the type of the PDE in (2). It is degen- 
erate at x = 0 and singular at x = +oo. It can be seen that r 2 := +oo is a natural boundary in 
any case, while the nature of r\ := 0, where the diffusion coefficient vanishes, depends on the 
parameters, a,b,c, see Appendix D. 

As Eq. (62) is in the adjoint form [cf. (2)], when c < 0, that is when one of the boundary 
points is an exit point and the other is natural (Appendix D), Feller theory guarantees that there 
is a unique solution without imposing any BC. In the other two cases discussed in Appendix D, 
instead, there is no uniqueness unless a BC at r, = 0 is prescribed; the boundary r 2 = +oo, being 
in all cases natural, does not require any BC to ensure uniqueness. 

In general, removing the potential term in (2) requires solving (usually by numerical means) 
a related ordinary differential equation. Setting v := u e~ hl , however, we are led again to a PDE 
like (1), 


u, = axu xx + (2a — c — bx) u x , 0 < x < +oo, t > 0. (63) 

Its Crank-Nicolson-type discretization is 


>(■ 


s i (u n + 1 - it')) 


J\ Lt j 


= MmW ~ O ~ “ “?!) + /WA. - »5) - («4) 


where B.,\ :=a.,u.,i, and 

r J ± 2 J ± 2 J ± 2 


six) := 


^ exp l/ 


b(r) 

a(r) 


dr\ — - x 


1 — - -- X 

a e a . 


(65) 


In Fig. 10, the evolution (growth) of a population governed by Eq. (62), with parameters 
a = 6, b = 2, and c = — 8 and an initial profile v(x, 0) = x exp {— x/5}, is shown. As c < a, 
the point x = 0 is an exit boundary, hence by Feller theory there is a unique solution without 
any BCs. We implemented the scheme described in Section II, using incomplete Choleski fac- 
torization to precondition it, and BCs of vanishing flux, which correspond to the true behavior of 
the aforementioned solution. In Fig. 1 1, we show the solution corresponding to the initial profile 
v(x, 0) = x 2 exp {— x 2 /5}. 

When 0 < c < a, the boundary point x = 0 is regular, and hence a BC is required there to 
recover a unique solution. If we seek, for instance, a solution with the boundary data u n Q = 1 (for 
all n), then we would implement 


si(u" +l - u '{ ) = ti[/3y 2 {u n 2 +1 - w" +1 ) - P l/2 (u1 +1 


i)+M«2-«;)-M“7-i)]’ 

(66) 
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Evolution of a population 



FIG. 10. Population evolution with initial profile u(.r,0) = x exp {— x/5) and parameters a = 6, b = 2, 
and c = —8. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 


that is 


u " +l U'l + ft(/?3/2 + ^1/2)] — M-^3/2 m 2 +1 — u l I 1 — M-(^3/2 + A/ 2 )] + 2/X/?i/2. (67) 

Given that fi ] / 2 is evaluated at the left end x = 0, the known term here vanishes, and the scheme 
coincides with that valid for the case with both BCs of vanishing flux. 


Evolution of a population 



FIG. 1 1 . Population evolution with initial profile u(x, 0) = x 2 exp { — jc 2 /5} and parameters a = 6, b = 2, 
and c = —8. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 


Numerical Methods for Partial Differential Equations DOI 10.1002/num 


NUMERICAL TREATMENT OF DEGENERATE DIFFUSION EQUATIONS 


827 


IV. SUMMARY 

Linear one-dimensional parabolic equations are considered in the pathological case that their dif- 
fusion coefficient vanishes on the boundary. The issue of uniqueness should be first addressed to 
decide whether BCs should be imposed, and, if so, which ones. Feller classification of the bound- 
ary points is used for this purpose. We solved problems having a unique solution without imposing 
arbitrary BCs. BCs which by Feller theory are known to be properties of such unique solutions are 
however used in a suitable algorithm. This is a preconditioned implicit finite-difference method 
on a suitable grid. The main properties of the algorithm are second-order accuracy, unconditional 
stability, and unconditional convergence to a constant as the time-step index tends to infinity. 

Problems as above but lacking uniqueness, generally speaking, require an interpretation, that 
is, should be better specified. A more precise or complete formulation should provide suitable 
additional data that make the full problem well posed, see Ref. [31], for instance. 

A recent analytical work [16] has been conducted on existence and uniqueness of solutions 
to degenerate elliptic and parabolic equations on open bounded connected domains in R" . The 
authors discuss when BCs have to be imposed, but not which ones they should or could be. 

APPENDIX A: FELLER CLASSIFICATION OF THE BOUNDARY POINTS 

The following results are due to Feller [13], but they are not always stated in a fully clear way in 
the literature. As this work is strongly based on Feller classification of the boundary points for 
the PDEs in (1) and (2), below we review and illustrate them. Defining the function 



where x 0 € (r \ , r 2 ) is arbitrary, which represents the Wronskian of the ordinary differential equa- 
tion a(x)y" + b(x)y' — 0 associated with the PDE in (1), up to the multiplicative constant W (jt 0 ), 
the boundary points were classified by Feller as follows [13]: 

• the boundary point rj, j = 1 or 2, is regular if 


a { (x)W l (x) <£ L l ((x 0 ,rj)) and W (x) a l (r)W 1 (r) dr e L l ((x 0 ,r y )); (A3) 


W(x) e L'((x 0 ,rj )) and a { (x)W 1 (x) e L l ((x 0 ,rj))’, 
• the boundary point rj is an exit point if 


(A2) 



• the boundary point rj is an entrance point if 



• the boundary point rj is natural in all other cases. 


In this article, we use extensively the function 



(A5) 
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and the flux related to Eq. (1), 

f(x, t ) := s(x) a{x) u x (x, t) = W _1 (x) u x (x, t). (A6) 

The classification of the boundary points allows one to determine which “lateral conditions” 
can be imposed on the Kolmogorov equations (1), (2), if any. Feller has established when a 
unique solution does exist, and when lateral conditions are needed to obtain it. Following Feller 
classification, we separate the following cases: 

a. Two natural boundaries. In this case, the initial-value problems for both PDEs, (1) and (2), 
have a unique solution, in the function space C°([ri , r 2 1) and L 1 ((r t , r 2 )), respectively, with- 
out imposing any lateral condition. In addition, Feller proved that in these cases the solution 
has zero flux at the boundaries, and such conditions can be explicitly prescribed (in a given 
algorithm), as they are indeed satisfied by the unique solution. 

b. An exit boundary and a natural boundary. Assume, for definiteness, that r\ is a natural 
boundary and r 2 is an exit boundary. Then, the initial- value problem associated with Eq. (1) 
has infinitely many solutions, whereas that for Eq. (2) has a unique solution. Such a unique 
solution has the property that a W v — >• 0 as r - > rf, and there is no need to impose any 
lateral condition here. This is instead needed to ensure uniqueness for Eq. (1). The solutions 
to (1) are characterized by the BC u(r 2 , t) = 0 (absorbing boundary). 

c. An entrance boundary and a natural boundary. For definiteness, let r, be a natural boundary 
and r 2 an entrance boundary. Then, the initial-value problem associated with Eq. (2) has 
infinitely many solutions, whereas that for Eq. (1) has a unique solution. The unique solution 
to (1) has the property that W _1 u x — > 0 as x — * rf and can be obtained without imposing 
any lateral condition. The latter is instead required to have a unique solution to (2). The 
solutions to (2) are characterized by the BC 

lim [( a(x ) v(x, t)) x — b(x) v(x, t)] = 0 

x ~* r 2 

(reflecting boundary). 

d. Two regular boundaries. When one or both the boundary points are regular, the initial-value 
problem for both equations, (1) and (2), has infinitely many solutions. Prescribing BCs is 
required for a solution of each of them to be unique. The solutions of Eq. (1) satisfy BCs of 
the form 


qj lim u(x,t) + (— 1 ) j Pj lim W 1 (x)u x (x,t ) = 0, 
whereas those of (2) satisfy 

qj lim W(x) a(x) v(x, t) + (— l) j pj lim {(a(x) v(x, t)) x — b(x) v(x, ?)} = 0, 

X— >Tj X ~^ r j 

for some constants pj, qj, j = 1,2. 

e. Any other case. When no boundary point is natural, two BCs need to be imposed to ensure 
uniqueness. 

A rather general property that is verified often is the following: when both boundaries 
are natural, or one is an entrance point and the other is natural, the initial-value problem for 
(1) has a unique solution. The boundary points being reflecting boundaries, the condition 
of vanishing flux can be used. In this case, the aforementioned unique solution attains a 
constant profile in the limit as t — +oo. 
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APPENDIX B: FINANCIAL MATHEMATICS 

Evaluating the Wronskian IV (S) = exp{— dxj, S 0 e (r ]5 r 2 ) = (0, +oo), where 

a(5) := S 2 and /3(S) := b — aS, we obtain 

I [ s b — ax ] b 

W (5) = exp | — y —dx^=kS a es, (Bl) 


and then: 

• fZ W(S)dS = f s j ) S"exp{l}dS. 

Setting j = 1, n = 0, we see that W(5) <£ L *(((), So)), hence t q = 0 is not a regular 
point. Setting j = 2, r 2 = +oo, we see that W(S) e L l ((S 0 , +oo)) if and only if a < — 1; 
for a > 0, r 2 — +oo is not a regular point. 

• /J a-'(S) W- l (S)dS = /'j S- 2 S- fl exp{-|}JS. 

Settingri = 0, we have a” 1 (5) W _1 (5) e L ! ((0, S 0 )); it follows that ri = 0 is not an exit 
point, but it satisfies the first condition to be an entrance point, see Section 2. Set r 2 = +oo. 
Then, a~'(5) W _1 (5) e L 1 ((S 0 , +oo)) if and only if a > —1; for a > 0, it follows that 
r 2 = +oo is not an exit point, but the first property for it to be an entrance point is satisfied. 


f 

Js 0 


-i 


(5) W 


-l 


( S ) f W(x)dx 
Js 0 


dS 


=*r 

JSn 


/sh* 0ex p{£}) rf * 


So S fl+2 exp { | } 

n n (S) 

Js 0 


dS 


D(S) 


dS. 


As D(5) — > +oo when S — 0 + and S — > +oo, L’Hopital’s rule can be applied to see that 
N(S)/D(S) ~ N'(S)/D'(S) = l /[ (a + 2)5 - b\, as 5 -> 0+ or 5 -* +oo, and hence 
N(S)/D(S) L'((5o, +oo)) and e L '((0,50)). This means that r\ = 0 is an entrance point, 

while r 2 = +oo is natural. 


APPENDIX C: WAVE PROPAGATION 


In this Appendix, we refer to Eqs. (56) and (57) for the lossy case and to Eq. (58) for the loss- 
less case. We apply Feller theory to classify the boundary points (in both coordinates, p and x). 
Let us start with Eq. (56). Defining, as usual, W(p) := exp{— ^ ds}, p 0 € (0, 1), where 
a(p) := j(l - p 2 ) 1 and b(p) := ±(l - p 2 ) 2 - \ a p, we have 


W(p) = exp 



2a s 

(1 - 5 2 ) 2 . 


ds 


(Cl) 


and setting u = s 2 , 


W(p) = exp 


-logs + 


a 


1 — s 2 


= k p 1 exp 


-‘P0 


a 

1 — p 2 


(C2) 
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Then, 


* /po W(p)dp = f' p J o p 1 expl^} dp. 

For r, = 0, W (p) turns out to be not integrable on (0, p 0 ); hence r j = 0 is not regular. 
For r 2 = 1, again W (p) £ L'((p 0 , 1)) and hence r 2 = 1 is not regular. 

* In a ”‘(P) W~\p)dp = f p J Q ( 1 - p 2 )~ 2 p expl --^} dp = [exp{— a-j-^)]p 0 . 

For r\ = 0, a~ l (p) W~ l {p) e L 1 ((0, p 0 )), hence r x = 0 is not an exit point but satisfies 
the first condition for being an entrance point. For r 2 = 1, a _1 (p) W 1 ( p ) e L'((p 0 , 1)), 
hence r 2 = 1 is not an exit point, but satisfies the first property to be an entrance point. 

* f p la-\p)W-\p)j p pQ W{s)dpds 

= fp 0 4 (! - P 2 )~ 2 P ex P {-7^2! fp 0 s ~ l ex P {7^2 }dp ds. 

Upon integrating by parts, we obtain [ e 1_ p 2 f p s 'e ] -P ds] p j 0 — [log p]f Q . 

Now, for r x = 0 both terms inside the squared parentheses diverge and hence a~ l (p)x 
x W~ l (p) f p W ( s ) ds is not integrable in (0, p 0 ), and the boundary r x — 0 is not an entrance 
point; hence, it is a natural boundary. For r 2 = 1 , the logarithm term is finite. To assess the 
behavior of the first term, we may use F’Hopital rule: 

a 

a rp _a f p s~ 1 ei -^ 2 ds 

lim e 1 ~p 2 / s~ 1 e ^ 1 ds = lim — 5 

J p Q P-*- 1“ e i- p 2 


It follows that a 1 
point. 


,. (1 -P 2 ) 2 A 

= lim — = 0. 

p^i- 2 p 2 

(p) W -1 (p) f p W(s)ds e L'CCpo, 1)) and hence r 2 = 1 is an entrance 


It is easy to see that the nature of the boundary points x = 0 and x = \ for the transformed 
PDE in (57), where a(x ) := x(l — x) 2 and b(x) := (1 — x) 2 — ax, is the same. We leave the 
check to the reader. 

Turning our attention to the PDE in (58), referring to a problem of wave propagation in a 
lossless random medium, we have 

, \ r a-* 2 ) 2 , 1 po 

W (p) = exp < — / — — dx\ = — , (C3) 

l J n x{ l-x-) 2 ) p 

and thus, up to a constant, 

* fpo W (p) d P = logPlpo- 

For /-| = 0, W(p) ^ L 1 (((), po)) and consequently r x = 0 is not regular. As for r 2 — 1, 
instead, W (p) e L l ((p Q , 1)); r 2 will be regular if the other condition will be verified: 

* Q ^'(p) d P = fp 0 77 ^ P d P = [7^72 ]«,- 

For r 1 = 0, the integral is finite, hence r, = 0 is not an exit point but it satisfies the 
first condition for being an entrance point. For r 2 = 1 , the integral diverges, hence r 2 = 1 
is neither a regular nor an entrance point, but the first condition for being an exit point is 
fulfilled. 

* W(p) f p a-\s)W~\s)dsdp = f 1 1 f p -^dsdp 

J P0 yr/ JpQ V V r J P0 P J P0 (1-S 2 F r 

= fp 0 f [7Z772 - 7 ^dp = k log p|‘ 0 + fl 7777^57 dp 
= k logp|‘ 0 + \_\ log 7^4- 
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This diverges at p = 1, hence W (p) a 1 
not an exit point and thus it is natural. 

< i lo S (s) ~ lo s (») - °- hence 

then r 2 = 0 is an entrance point. 


(s) W x (s)ds £ L ] ((p 0 , 1)) and r 2 = 1 is 

f p - ds dp = log ( rfp, and 

Jp 0 i ^ (l_ p 2)2 b \poJ r 

( P ) W~\p)f p Q W(s)ds e L'((0,p 0 )), and 


APPENDIX D: GENETICS 

The classification of the boundary points yield the properties of solutions as reported by Feller in 
Ref. [3]. We shall see that the boundary x = +oo is natural and hence no BC has to be imposed 
there. The boundary x = 0 is of a different kind according to whether c < 0, 0 < c < a, or 
c > a. In Ref. [3] Feller showed that 

i. if c < 0, the initial value u(x, 0) can be prescribed arbitrarily and it determines uniquely a 
solution without imposing any BC. Such solution preserves positivity and has a decreasing 
norm; 

ii. if 0 < c < a, the case is similar to that of a regular equation, in that there exists a solution 
to the initial value problem, which is positive, preserves its norm, and is defined by the BC 
of vanishing flux at x = 0. There are infinitely many other solutions that preserve positivity 
and possess a decreasing norm, among which only one satisfies the condition w(0, t) < oo. 
Moreover, there exists a unique solution with initial values identically equal to zero and a 
flux prescribed at x = 0; this solution is integrable but in general not bounded; 

iii. if c > a, there exists a solution to the initial- value problem, which is positive and norm 
preserving, and vanishing along with its flux at x = 0. This means that x = 0 is both 
an absorbing and a reflecting barrier, and thus no homogeneous BC can be imposed there. 
Solutions with an arbitrarily prescribed flux at x = 0 do exist, which may take even negative 
values or have an increasing norm, for some values of t. 

We can recover some of Feller’s findings studying directly the nature of the boundary points. 
Correspondingly to Eq. (62), we have 


W(x) = exp 



= k x 


(Dl) 


with x 0 e (0, Too) arbitrary, and a ( 5 ) := as, /3(s) := bs T c. Then, 

• fxQ W(x)dx = f'i x~^e~a x dx. 

For r, = 0, W(jc) e L'((0, x 0 )) if and only if c < a. For r 2 = Too, W(x) e 
L x ((jt 0 . Too)) if b > 0. 

• f^ J a -1 (x) W -1 (x) dx = j' x ^ (ax)~ ] x 'a ea x dx. 

For r\ = 0, a~ l (x) W~ l (x) e L l ((0,x 0 )) if and only if c - >0, i.e., c > 0; it follows 
that r j = 0 is regular if and only if 0 < c < a\ when c < 0, instead, the first condition 
for the boundary r x = 0 to be an exit point is satisfied. For r 2 = Too, then if b > 0, 
a _1 (x) W~ x (x) £ L l ((jt 0 , Too)); hence, r 2 = Too is not regular. 

• f^ j W(x) o' 1 (^-) W _1 (5) ds dx = f' j x exp {— | x} f* (as)~ l s« exp s} ds dx. 
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To check the second condition upon which a given boundary is an exit point, we integrate 
repeatedly by parts in f* (as)-' s '<• e « ' ds, integrating with respect to x and differentiating 

the exponential, j'J x~«e~^ x x^e^ x dx = r ; — x 0 . 

For r i = 0, W(x) f x cr -1 ^) e L‘((0,xo)), hence r\ = 0 is an exit point 

when c < 0. For r 2 = +oc, instead, this function is not integrable in (x 0 , +oo) and hence 
r 2 = +oo is not an exit point. 

• fj a _1 (x)W _1 (x) f* W (s) ds dx = (ax)~ 1 x“ e% x f* s~^ e ~a s ds dx. 

Integrating again by parts, we obtain x°~ 1 ea x x~^ +1 ^e~a s dx = r, — x 0 . 

For r\ = 0, cr 1 (x) W _1 (x) f x W(s)ds e L l ((0, x 0 )); as the first condition upon which 
the boundary is an entrance point was satisfied for c > 0, it follows that r, = 0 is an entrance 
point for c < 0. As for 0 < c < a the boundary is regular, = 0 is an entrance point for 
c > a. For r 2 = +oo, instead, this function is not integrable in (x 0 , +oo), hence r 2 = +oo 
is not an entrance point and thus it is a natural boundary in any case. 
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